%% process_pde_NS_t
function pde = process_pde_NS_t(index)
switch index
    case 1
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y,t) (x.^2.*y.^2 + exp(-y)).*cos(2*pi*t);
        pde.u2 = @(x,y,t) (-(2/3)*x.*y.^3 + 2 - pi*sin(pi*x)).*cos(2*pi*t);
        pde.p = @(x,y,t) -(2-pi*sin(pi*x)).*cos(2*pi*y).*cos(2*pi*t);
        pde.u1dx = @(x,y,t) 2*x.*y.^2.*cos(2*pi*t);
        pde.u1dy = @(x,y,t) -cos(2*pi*t).*(-2*y.*x.^2 + exp(-y));
        pde.u2dx = @(x,y,t) -cos(2*pi*t).*((2*y.^3)./3 + pi^2*cos(x*pi));
        pde.u2dy = @(x,y,t) -2*x.*y.^2.*cos(2*pi*t);
        pde.pdx = @(x,y,t) pi^2*cos(2*pi*t).*cos(pi*x).*cos(2*pi*y);
        pde.pdy = @(x,y,t) -2*pi*cos(2*pi*t).*sin(2*pi*y).*(pi*sin(pi*x) - 2);
        pde.f1 = @(x,y,t) cos(2*pi*t).^2.*(-2*y.*x.^2 + exp(-y)).*((2*x.*y.^3)./3 + pi*sin(pi*x) - 2) - nu*(cos(2*pi*t).*(2*x.^2 + exp(-y)) + 2*y.^2.*cos(2*pi*t)) - 2*pi*sin(2*pi*t).*(exp(-y) + x.^2.*y.^2) + 2*x.*y.^2.*cos(2*pi*t).^2.*(exp(-y) + x.^2.*y.^2) + pi^2*cos(2*pi*t).*cos(pi*x).*cos(2*pi*y);
        pde.f2 = @(x,y,t) 2*pi*sin(2*pi*t).*((2*x.*y.^3)./3 + pi*sin(pi*x) - 2) - cos(2*pi*t).^2.*((2*y.^3)./3 + pi^2*cos(x*pi)).*(exp(-y) + x.^2.*y.^2) - nu*(pi^3*cos(2*pi*t).*sin(pi*x) - 4*x.*y.*cos(2*pi*t)) - 2*pi*cos(2*pi*t).*sin(2*pi*y).*(pi*sin(pi*x) - 2) + 2*x.*y.^2.*cos(2*pi*t).^2.*((2*x.*y.^3)./3 + pi*sin(pi*x) - 2);
        pde.domain = [0, 1, -0.25, 0, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    case 2
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y,t) pi*sin(2*pi*y).*sin(pi*x).^2.*sin(t);
        pde.u2 = @(x,y,t) -pi*sin(2*pi*x).*sin(pi*y).^2.*sin(t);
        pde.p = @(x,y,t) 10*cos(pi*x).*sin(pi*y).*sin(t);
        pde.u1dx = @(x,y,t) 2*pi^2*cos(pi*x).*sin(pi*x).*sin(2*pi*y).*sin(t);
        pde.u1dy = @(x,y,t) 2*pi^2*cos(2*pi*y).*sin(pi*x).^2.*sin(t);
        pde.u2dx = @(x,y,t) -2*pi^2*cos(2*pi*x).*sin(pi*y).^2.*sin(t);
        pde.u2dy = @(x,y,t) -2*pi^2*cos(pi*y).*sin(2*pi*x).*sin(pi*y).*sin(t);
        pde.pdx = @(x,y,t) -10*pi*sin(pi*x).*sin(pi*y).*sin(t);
        pde.pdy = @(x,y,t) 10*pi*cos(pi*x).*cos(pi*y).*sin(t);
        pde.f1 = @(x,y,t) pi*sin(pi*x).^2.*sin(2*pi*y).*cos(t) - nu*(2*pi^3*cos(pi*x).^2.*sin(2*pi*y).*sin(t) - 6*pi^3*sin(pi*x).^2.*sin(2*pi*y).*sin(t)) - 10*pi*sin(pi*x).*sin(pi*y).*sin(t) + 2*pi^3*cos(pi*x).*sin(pi*x).^3.*sin(2*pi*y).^2.*sin(t).^2 - 2*pi^3*cos(2*pi*y).*sin(pi*x).^2.*sin(2*pi*x).*sin(pi*y).^2.*sin(t).^2;
        pde.f2 = @(x,y,t) nu*(2*pi^3*cos(pi*y).^2.*sin(2*pi*x).*sin(t) - 6*pi^3*sin(2*pi*x).*sin(pi*y).^2.*sin(t)) - pi*sin(2*pi*x).*sin(pi*y).^2.*cos(t) + 10*pi*cos(pi*x).*cos(pi*y).*sin(t) + 2*pi^3*cos(pi*y).*sin(2*pi*x).^2.*sin(pi*y).^3.*sin(t).^2 - 2*pi^3*cos(2*pi*x).*sin(pi*x).^2.*sin(pi*y).^2.*sin(2*pi*y).*sin(t).^2;
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    case 3
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y,t) 10.*x.^2.*y.*cos(t).*(2.*y - 1).*(x - 1).^2.*(y - 1);
        pde.u2 = @(x,y,t) -10.*x.*y.^2.*cos(t).*(2.*x - 1).*(x - 1).*(y - 1).^2;
        pde.p = @(x,y,t) cos(t).*(20.*x - 10).*(2.*y - 1);
        pde.u1dx = @(x,y,t) 20*x.*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2*x - 2).*(2*y - 1).*(y - 1);
        pde.u1dy = @(x,y,t) 10*x.^2.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 20*x.^2.*y.*cos(t).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2*y - 1).*(x - 1).^2;
        pde.u2dx = @(x,y,t) -10*y.^2.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 - 20*x.*y.^2.*cos(t).*(x - 1).*(y - 1).^2 - 10*x.*y.^2.*cos(t).*(2*x - 1).*(y - 1).^2;
        pde.u2dy = @(x,y,t) -20.*x.*y.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 - 10*x.*y.^2.*cos(t).*(2*x - 1).*(2*y - 2).*(x - 1);
        pde.pdx = @(x,y,t) 20.*cos(t).*(2*y - 1);
        pde.pdy = @(x,y,t) 2*cos(t).*(20*x - 10);
        pde.f1 = @(x,y,t) 20*cos(t).*(2*y - 1) - nu*(40*x.^2.*cos(t).*(x - 1).^2.*(y - 1) + 20*x.^2.*cos(t).*(2*y - 1).*(x - 1).^2 + 40*x.^2.*y.*cos(t).*(x - 1).^2 + 20*x.^2.*y.*cos(t).*(2*y - 1).*(y - 1) + 20*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 40*x.*y.*cos(t).*(2*x - 2).*(2*y - 1).*(y - 1)) - 10*x.^2.*y.*sin(t).*(2*y - 1).*(x - 1).^2.*(y - 1) - 10*x.*y.^2.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2.*(10*x.^2.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 20*x.^2.*y.*cos(t).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2.*y - 1).*(x - 1).^2) + 10*x.^2.*y.*cos(t).*(2*y - 1).*(20*x.*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1) + 10*x.^2.*y.*cos(t).*(2*x - 2).*(2*y - 1).*(y - 1)).*(x - 1).^2.*(y - 1);
        pde.f2 = @(x,y,t) 2*cos(t).*(20*x - 10) + nu*(40*y.^2.*cos(t).*(x - 1).*(y - 1).^2 + 20*y.^2.*cos(t).*(2*x - 1).*(y - 1).^2 + 40*x.*y.^2.*cos(t).*(y - 1).^2 + 20*x.*y.^2.*cos(t).*(2*x - 1).*(x - 1) + 20*x.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 + 40*x.*y.*cos(t).*(2*x - 1).*(2*y - 2).*(x - 1)) + 10*x.*y.^2.*sin(t).*(2*x - 1).*(x - 1).*(y - 1).^2 - 10*x.^2.*y.*cos(t).*(2*y - 1).*(x - 1).^2.*(y - 1).*(10*y.^2.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 + 20*x.*y.^2.*cos(t).*(x - 1).*(y - 1).^2 + 10*x.*y.^2.*cos(t).*(2*x - 1).*(y - 1).^2) + 10*x.*y.^2.*cos(t).*(2*x - 1).*(20*x.*y.*cos(t).*(2*x - 1).*(x - 1).*(y - 1).^2 + 10*x.*y.^2.*cos(t).*(2*x - 1).*(2*y - 2).*(x - 1)).*(x - 1).*(y - 1).^2;
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    case 4
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y,t) 5*(x.^2-2*x.^3+x.^4).*(2*y-6*y.^2+4*y.^3).*cos(2*pi*t);
        pde.u2 = @(x,y,t) -5*(y.^2-2*y.^3+y.^4).*(2*x-6*x.^2+4*x.^3).*cos(2*pi*t);
        pde.p = @(x,y,t) 5*(2*x-1).*(2*y-1).*cos(2*pi*t);
        pde.u1dx = @(x,y,t) x.*y.*cos(t.*pi.*2.0).*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*2.0e+1;
        pde.u1dy = @(x,y,t) x.^2.*cos(t.*pi.*2.0).*(x-1.0).^2.*(y.*-1.2e+1+y.^2.*1.2e+1+2.0).*5.0;
        pde.u2dx = @(x,y,t) y.^2.*cos(t.*pi.*2.0).*(y-1.0).^2.*(x.*-1.2e+1+x.^2.*1.2e+1+2.0).*-5.0;
        pde.u2dy = @(x,y,t) x.*y.*cos(t.*pi.*2.0).*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*-2.0e+1;
        pde.pdx = @(x,y,t) cos(t.*pi.*2.0).*(y.*2.0-1.0).*1.0e+1;
        pde.pdy = @(x,y,t) cos(t.*pi.*2.0).*(x.*1.0e+1-5.0).*2.0;
        pde.f1 = @(x,y,t) cos(t.*pi.*2.0).*(y.*2.0e+1-1.0e+1)-nu.*cos(t.*pi.*2.0).*(y.*2.0-1.0).*(-y+x.^2.*y.^2.*6.0+x.*y.*6.0-x.*y.^2.*6.0-x.^2.*y.*6.0+x.^2.*3.0-x.^3.*6.0+x.^4.*3.0+y.^2).*2.0e+1-x.^2.*y.*pi.*sin(t.*pi.*2.0).*(x-1.0).^2.*(y.*-3.0+y.^2.*2.0+1.0).*2.0e+1+x.^3.*y.^2.*cos(t.*pi.*2.0).^2.*(x-1.0).^2.*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).^2.*2.0e+2-x.^3.*y.^2.*cos(t.*pi.*2.0).^2.*(x-1.0).^2.*(y-1.0).^2.*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-1.2e+1+y.^2.*1.2e+1+2.0).*5.0e+1;
        pde.f2 = @(x,y,t) cos(t.*pi.*2.0).*(x.*2.0e+1-1.0e+1)+nu.*cos(t.*pi.*2.0).*(x.*2.0-1.0).*(-x+x.^2.*y.^2.*6.0+x.*y.*6.0-x.*y.^2.*6.0-x.^2.*y.*6.0+x.^2+y.^2.*3.0-y.^3.*6.0+y.^4.*3.0).*2.0e+1+x.*y.^2.*pi.*sin(t.*pi.*2.0).*(y-1.0).^2.*(x.*-3.0+x.^2.*2.0+1.0).*2.0e+1+x.^2.*y.^3.*cos(t.*pi.*2.0).^2.*(y-1.0).^2.*(x.*-3.0+x.^2.*2.0+1.0).^2.*(y.*-3.0+y.^2.*2.0+1.0).*2.0e+2-x.^2.*y.^3.*cos(t.*pi.*2.0).^2.*(x-1.0).^2.*(y-1.0).^2.*(x.*-1.2e+1+x.^2.*1.2e+1+2.0).*(y.*-3.0+y.^2.*2.0+1.0).*5.0e+1;
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    case 5
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y,t) (x.^2-2*x.^3+x.^4).*(2*y-6*y.^2+4*y.^3).*cos(t);
        pde.u2 = @(x,y,t) -(2*x-6*x.^2+4*x.^3).*(y.^2-2*y.^3+y.^4).*cos(t);
        pde.p = @(x,y,t) (x.^2+y.^2-2/3).*cos(t);
        pde.u1dx = @(x,y,t) x.*y.*cos(t).*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*4.0;
        pde.u1dy = @(x,y,t) x.^2.*cos(t).*(x-1.0).^2.*(y.*-1.2e+1+y.^2.*1.2e+1+2.0);
        pde.u2dx = @(x,y,t) -y.^2.*cos(t).*(y-1.0).^2.*(x.*-1.2e+1+x.^2.*1.2e+1+2.0);
        pde.u2dy = @(x,y,t) x.*y.*cos(t).*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*-4.0;
        pde.pdx = @(x,y,t) x.*cos(t).*2.0;
        pde.pdy = @(x,y,t) y.*cos(t).*2.0;
        pde.f1 = @(x,y,t) x.*cos(t).*2.0-nu.*cos(t).*(y.*2.0-1.0).*(-y+x.^2.*y.^2.*6.0+x.*y.*6.0-x.*y.^2.*6.0-x.^2.*y.*6.0+x.^2.*3.0-x.^3.*6.0+x.^4.*3.0+y.^2).*4.0-x.^2.*y.*sin(t).*(x-1.0).^2.*(y.*-3.0+y.^2.*2.0+1.0).*2.0+x.^3.*y.^2.*cos(t).^2.*(x-1.0).^2.*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).^2.*8.0-x.^3.*y.^2.*cos(t).^2.*(x-1.0).^2.*(y-1.0).^2.*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-1.2e+1+y.^2.*1.2e+1+2.0).*2.0;
        pde.f2 = @(x,y,t) y.*cos(t).*2.0+nu.*cos(t).*(x.*2.0-1.0).*(-x+x.^2.*y.^2.*6.0+x.*y.*6.0-x.*y.^2.*6.0-x.^2.*y.*6.0+x.^2+y.^2.*3.0-y.^3.*6.0+y.^4.*3.0).*4.0+x.*y.^2.*sin(t).*(y-1.0).^2.*(x.*-3.0+x.^2.*2.0+1.0).*2.0+x.^2.*y.^3.*cos(t).^2.*(y-1.0).^2.*(x.*-3.0+x.^2.*2.0+1.0).^2.*(y.*-3.0+y.^2.*2.0+1.0).*8.0-x.^2.*y.^3.*cos(t).^2.*(x-1.0).^2.*(y-1.0).^2.*(x.*-1.2e+1+x.^2.*1.2e+1+2.0).*(y.*-3.0+y.^2.*2.0+1.0).*2.0;
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    case 6
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y,t)-sin(pi*x).^2.*sin(2*pi*y).*cos(t);
        pde.u1dx = @(x,y,t)pi.*cos(x.*pi).*sin(x.*pi).*sin(y.*pi.*2.0).*cos(t).*-2.0;
        pde.u1dy = @(x,y,t)pi.*cos(y.*pi.*2.0).*sin(x.*pi).^2.*cos(t).*-2.0;
        pde.u2 = @(x,y,t)sin(2*pi*x).*sin(pi*y).^2.*cos(t);
        pde.u2dx = @(x,y,t)pi.*cos(x.*pi.*2.0).*sin(y.*pi).^2.*cos(t).*2.0;
        pde.u2dy = @(x,y,t)pi.*cos(y.*pi).*sin(x.*pi.*2.0).*sin(y.*pi).*cos(t).*2.0;
        pde.p = @(x,y,t)(x.*y-0.25).*cos(t);
        pde.pdx = @(x,y,t)y.*cos(t);
        pde.pdy = @(x,y,t)x.*cos(t);
        pde.f1 = @(x,y,t)nu.*(pi.^2.*cos(x.*pi).^2.*sin(y.*pi.*2.0).*cos(t).*2.0-pi.^2.*sin(x.*pi).^2.*sin(y.*pi.*2.0).*cos(t).*6.0)+y.*cos(t)+sin(x.*pi).^2.*sin(y.*pi.*2.0).*sin(t)+pi.*cos(x.*pi).*sin(x.*pi).^3.*sin(y.*pi.*2.0).^2.*cos(t).^2.*2.0-pi.*cos(y.*pi.*2.0).*sin(x.*pi).^2.*sin(x.*pi.*2.0).*sin(y.*pi).^2.*cos(t).^2.*2.0;
        pde.f2 = @(x,y,t)-nu.*(pi.^2.*cos(y.*pi).^2.*sin(x.*pi.*2.0).*cos(t).*2.0-pi.^2.*sin(x.*pi.*2.0).*sin(y.*pi).^2.*cos(t).*6.0)+x.*cos(t)-sin(x.*pi.*2.0).*sin(y.*pi).^2.*sin(t)+pi.*cos(y.*pi).*sin(x.*pi.*2.0).^2.*sin(y.*pi).^3.*cos(t).^2.*2.0-pi.*cos(x.*pi.*2.0).*sin(x.*pi).^2.*sin(y.*pi).^2.*sin(y.*pi.*2.0).*cos(t).^2.*2.0;
        pde.domain = [0 1 0 1 0 1];
        pde.bdary = [0 0 0 0];
    otherwise
        error("Invalid pde index.");
end
fprintf("pde index: %d\n", pde.index);
fprintf("pde domain: xmin,xmax,ymin,ymax,tmin,tmax\n");
disp(pde.domain);
fprintf("pde boundary: bottom,right,top,left (0=Dirichlet,1=Neumann)\n");
disp(pde.bdary);
if isfield(pde, "bdary_p")
    fprintf("pde boundary of p: (0=Dirichlet,1=Neumann)\n");
    disp(pde.bdary_p);
end
end